Modeling of the combined dynamics of leptospirosis transmission and seroconversion in herds

Leptospirosis is a zoonotic disease-causing illness in both humans and animals resulting in related economic impacts due to production loss as well as prevention and control efforts. Several mathematical models have been proposed to study the dynamics of infection but none of them has so far taken into account the dynamics of seroconversion. In this study, we have developed a general framework, based on the kinetic model for animal leptospirosis, that combines both the antibody (exposure marker) and infection dynamics to simultaneously follows both seroconversion and infection status of leptospirosis in a herd population. It is a stochastic compartmental model (for transition rates) with time delay (for seroconversion) which describes the progression of infection by a SEIRS (susceptible, exposed, infected, removed and susceptible) approach and seroconversion by four-state antibody kinetics (antibody negative and three antibody positive states of different antibody levels). The model shows that it is possible to assess and follow both seroconversion and infection status through the prism of diagnostic testing. Such an approach of combined kinetics could prove very useful to assist the competent authorities in their analyzes of epidemic situations and in the implementation of strategies for controlling and managing the associated risks.

examination as it produces high sensitivity and specificity results 7 . In addition, MAT is used to identify the circulated Leptospira serogroups 8-10 by simultaneously detecting both immunoglobulin M (IgM) and immunoglobulin G (IgG), which are classes of agglutinating antibodies 11,12 . Nonetheless, ELISA can differentiate IgM or IgG antibodies 10,13 . To detect the antigen, PCR is an alternative. This method can directly detect leptospiral DNA in the samples. However, it cannot identify the etiological serovars 12 . The isolation and identification of Leptospira is possible with the culture method. Nevertheless, the method is difficult for many reasons, for examples, type of samples and the timing of taken samples 12 . Therefore, all currently available leptospirosis infection diagnostic tests could not provide a definite indication.
To study the infection dynamics of leptospirosis, several mathematical models have been proposed. A Susceptible-Infectious (SI) epidemiological model was used for the spread of leptospirosis in livestock by varying periodic parameters 14 and for the leptospire dynamics and control in the Norway rat 15 . A Susceptible-Infectious-Retired (SIR) model was used to study the transmission of leptospirosis between human and animal populations 16,17 . This SIR model was then used to understand the epidemiology of leptospirosis in cattle 18 . To improve the description of the dynamics of infection with Leptospira, a SEIR (SIR plus the "Exposed" class) model was used to study the optimal control of disease outbreak 19 . However, those models mainly focused on the sensitivity of the transmission parameters that are related to the number of infected animals based on the simulation results. Nevertheless, the time series of infected animals are not well defined due to the limitations of current diagnosis methods. To follow the spread of the disease in the animal population, serological diagnostics are widely used and accepted. However, to our knowledge, the seroconversion kinetics of leptospirosis have never been taken into account in previous mathematical models cited above. Indeed, by supplementing those models on the dynamics of the infection by tracking the level of antibodies (exposure marker or probe) makes it possible to have a more detailed view of the state of the circulation of pathogens in the considered population and to improve the interpretations of diagnostic tests.
In this study, we aim to assess the epidemiological correlations between the seroconversion dynamics, as obtained from serological diagnostic tests, and infectious status of a population of animals during an outbreak. To this end, we built a general framework based on the kinetic model for animal leptospirosis that combines both the infection and the antibody dynamics. Using a compartmental model with time delay, the model simultaneously follows both infection and seroconversion dynamics. The states of the infection and antibody classes of the population are described as a function of the basic reproduction number allowing the correlation between the number of infections (epizootic size) and the prevalence of antibody-positive individuals.

Results
Seroconversion dynamical model. A first key result of this study is the description of the general framework of the seroconversion dynamical model consisted of a combination of the infection and antibody dynamics. The kinetic schemes of the infection dynamics, antibody dynamics and combined model are shown in Fig. 1 and the equations are described in the "Methods" section.
The infection dynamics was explained following the compartmental Susceptible-Exposed-Infectious-Recovered-Susceptible (SEIRS) model. Rather than a SIR model like in some works cited above, we have chosen to use a SEIRS model to allow a stage (compartment E) infected but not yet infectious and which may or may not already be carrying antibodies. Indeed, the development and appearance of antibodies taking place during the incubation phase, such a stage is quite plausible and possible. The SEIRS model is richer and can be reduced to a SIRS model when the duration of stage E turns out to be very short; the reverse not being possible. Moreover, as we used a stochastic SEIRS model with residence times in the stages obtained from the constrained distributions of Eq. (6), epidemiological situations with or without stage E are statistically possible.
Infectious animals can transmit the disease to susceptible animals that progress to exposed animals with a transmission rate . Exposed animals became infectious after an incubation period of mean duration t e = 1/ν . The mean duration of infection is t i = 1/α . Infectious animals progress to recovered and immune state, which became susceptible again after a mean duration of immunity of t r = 1/γ . All animals die at a mortality rate µ considered as the population renewal rate. The set of equations related to the dynamics of infection were provided in the Supplement Information (S1. Methods).
Regarding the antibody kinetics, the model is based on the antibody level of animals (mainly cattle here). Primarily exposed animals developed antibody levels from negative to positive in both IgM and IgG antibodies and later only positive IgG antibodies. IgM antibodies stay for 3-5 weeks (in cattle) 20 while IgG antibodies, appearing at the same time or just after IgM antibodies, stay for a much longer time 20 . Subsequently, the animal IgG antibodies decrease at a rate ω 1 and became negative again in both IgM and IgG antibodies. Afterwards, re-infection of those animals rose their IgG antibodies to a higher level 21,22 (IgG + state).
The infection and the antibody kinetics are combined to study the seroconversion dynamics. When susceptible animals were first infected with leptospires, they enter exposed state with negative antibodies. Subsequently, exposed animals develop IgM and IgG antibodies at a rate ν 0 , and then progress to the infectious state at a rate ν ea . However, infectious animals that recovered still carried IgG antibodies during the immunity period and became susceptible again with negative results in both IgM and IgG antibodies at a rate γ.
Model outcomes. The average number of infectious animals (i.e., IA 2 , IA 1 and IA 3 compartments) are shown in Fig. 2 for the basic reproduction number R 0 equal to 1.5, 2.5, and 5.0. All results of the 10 compartments are shown in Fig. S1. As a check, the comparisons between the combined model and the reduced SEIRS model are shown in Fig. S2. At the beginning, there was no circulation of leptospirosis infection, the epidemic curve begins to increase and then oscillates around a plateau. The results showed that the epidemic curve could be divided into two phases: a growth phase and a stationary state. During the growth phase, the total infectious www.nature.com/scientificreports/ animals is mainly composed of IA 2 and IA 1 . After that, the total infected population consists of the number of IA 3 , IA 2 , and IA 1 , where most of the infectious animals came from IA 3 . For R 0 > 1 , the first peak of infectious animals comes from the number of IA 2 , carrying IgM and IgG antibodies, followed by IA 1 and IA 3 carrying IgG and high-level of IgG antibodies, respectively. The time interval between first peaks of each infectious compartments changed with R 0 . For all value of R 0 , the time interval between peaks of IA 2 -IA 1 is about twice of the IgM duration ( 1/ω 2 ). The time interval between peaks of IA 2 -IA 3 and IA 1 -IA 3 decreased as R 0 increases. Increasing IgM duration leads to increase the time interval between peaks of IA 2 -IA 3 , and IA 1 -IA 3 , but does not affect the time interval between peaks of IA 2 -IA 1 (data not shown). The cross-correlations of the monthly number of IA 2 , IA 1 and IA 3 are plotted in Fig. S3. A positive strong correlation between IA 1 and IA 2 with a positive 4-month lag was found for all value of R 0 , excepted for R 0 = 1.5, the lag time was 3 months. This 4-month lag corresponds to the time interval between the peaks of the number of IA 2 and IA 1 (Fig. 2). Overall, the correlation between IA 3 -IA 2 and IA 3 -IA 1 were similar excepted for R 0 = 1.5 . From a − 1.5 to 1.5 years lag, the negative correlations between IA 3 -IA 2 and IA 3 -IA 1 were found, which implied that the high number of IA 2 and IA 1 related to the low number of IA 3 . The simulation results can be used to forecast the number of IA 1 and IA 3 at later times from the number of IA 2 in the past (occurring first); that is to say, when the number of new infectious IA 2 was observed, the number of IA 1 or IA 3 can be predicted. Prior to the stationary state, Fig. 2 shows that there is no general relation in the numbers of infected animals as a function of R 0 .

Characteristics of infectious.
Diagnosis tests such as ELISA can be used to inform the public health about the situation concerning the number of infected in each animal herd. The number of infected animals is derived and classified as the numbers of TR 1 , TR 2 , and TR 3 (as shown in Fig. 6 in "Methods" section). The total number of infectious and that of TR 1 , TR 2 , and TR 3 are shown in Fig. 3. Both Figs. 2 and 3 exhibit similar trends as a function of time. During the growth phase, the total number of infectious animals is mainly composed of the number of TR 1 and TR 2 , afterward it tends to be mainly TR 2 and TR 3 . The cross-correlations between the number of TR i are shown in Fig. S4. Overall, the patterns of the cross-correlations between the number of infectious animals were similar to the number of TR i . The positive strong correlations between TR 1 and TR 2 were found to decrease as R 0 increase. This indicates that the number of TR 2 can be forecasted from the number of TR 1 with a lag time L 1 . Figure 4 shows that L 1 tends to decrease a bit with R 0 . However, positive weak correlations between the number  Table 1. TR i stands for "Test Result" of a diagnostic test. Arrows tilted up correspond to the natural mortality rate μ (identical for all arrows but just indicated for some).

Scientific Reports
| (2022) 12:15620 | https://doi.org/10.1038/s41598-022-19833-x www.nature.com/scientificreports/ of TR 3 and TR 1 were observed at lag times ( L 2 ) of 50, 44, 37 and 35 months for R 0 = 1.5, 2.5, 5.0 and 7.5, respectively. As shown in Fig. 4, L 2 decreases with R 0 . The comparison of results between the number of infectious IA i and the number of TR i for R 0 = 1.5, 2.5, and 5.0 is shown in Fig. S5. The graph indicates that the number of TR i perfectly parallel the number of infectious animals (with, TR 1 ∼ IA 2 ) all the time both in the growth and stationary phases. Therefore, the number of TR i can be used to extract the proportion of infectious animals any time. To estimate a cut-off time between growth and stationary phases, the relaxation time (time to a stationary state or a cut-off time ( t c )) is defined as the time reached at the 1st minimum of TR i after the 1st peak of that curve. The average number of TR i at stationary state is calculated as the average number over the period since the relaxation time to the end of simulation time (30 years). Values of t c for TR i are provided in Table S2.
The average number of TR i at stationary state was used to calculate the probability p i of finding antibody positive animals in TR i as a function of R 0 (Fig. 5, left panel). Clearly, p i is equal to zero for R 0 ≤ 1 and increases with R 0 . However, the probability q i of finding infected individuals among the positive TR i showed no variation with R 0 (Fig. 5, right panel) as predicted from Eq. (10) for the stationary state expression of q i,s .

Discussion and conclusions
Leptospirosis is a zoonotic disease-causing illness in both humans and animals resulting in related economic impacts due to production loss as well as prevention and control efforts 23 . Several mathematical models have been constructed to study the dynamics of infection using compartment models, such as the SIR [16][17][18] and SEIR models 19 . These models did not consider the seroconversion dynamics. In this study, a general framework that combines both antigen and antibody dynamics was constructed to assess the state of leptospirosis in a population of animals using diagnostic tests. To our knowledge, this model is the first compartment model built to access leptospirosis progression profile in animals considering the IgM and IgG antibody kinetics and outcomes of diagnosis tests.
It is well known that the clinical signs in animals are mild or asymptomatic 24 . Thus, the situation of leptospirosis outbreak is only described through the test results from sampling population. In our framework, a novel approach to explore the dynamics of leptospirosis in animals was proposed with the number of infected animals derived from the test results.
Using IgM-ELISA, for example, with sensitivity 86% ( Se 1 ) and specificity 84% ( Sp 1 ) 25 , the fraction of positive T +,i makes it possible to calculate the prevalence p 1 or the true positive antibody fraction of animals. This is helpful to uncover the actual epidemic situation rather than relying on only the serological test results. Furthermore, the basic reproduction number ( R 0 ) is estimated using Eq. (9). In general, R 0 is a useful tool for describing and predicting the dynamics of infection. In mathematical models 15,26 , R 0 is defined from the transmission rate, a parameter that is difficult to access. In this study, R 0 can be prior calculated to infer the transmission rate. Subsequently, after the parameters in Eq. (10) are known, the fractions of infected states among the positive tests can be calculated. Thus, our model is useful to elucidate the picture of leptospirosis dynamics either from the parameters of Eq. (10) or from the data set from the field. Furthermore, the cross-correlation of the number of TR i was studied to provide better understanding on the time lag of infection. Animals with positive antibody is anticipated in later time after the outbreak. Using the time lag, for example, if the number of TR 1 is estimated, the number of TR 2 and TR 3 can be then predicted. Reversely, if the high titer level of IgG is known by other means the number of TR 3 can be calculated and the TR 1 and TR 2 in the past can be tracked back. The time lag output is based on the input parameters and is flexible to change. Thus, this model is also applicable and can be adapted for leptospirosis in humans once data collected from human outbreaks and relevant diagnostic tests is available.
The developed model that we have just described is based on only on few assumptions and the outcomes depend on the quality of the parameters used. The model considered a closed system with a homogeneously well-mixed population with the parameters used obtained from the literature (see Table 1). To take account for the uncertainties in the parameters we used distributions to sample the parameters. It was assumed that the parameters were constant over the entire study period. The model can be embellished in several directions including age structure and precise renewal pattern of populations, spatio-temporal aspects or other factors that may have an impact on the transmission of leptospires like weather conditions.
In conclusion, we developed and constructed a general and generic modeling framework allowing to describe the transmission dynamics of animal leptospirosis in which antibody kinetics is taken into account. The model also makes it possible to describe the outbreak situations through the prism of diagnostic tests. This proposed new approach could prove useful to the competent authorities in their analyzes of real outbreak situations and, consequently, in the implementation of disease control and associated risk management strategies.

Methods
Seroconversion dynamical model. The kinetics of antibody levels in an animal with leptospirosis can be described as follows (Fig. 1, left bottom). A 0 represents the numbers of naïve animals that have never encountered leptospires and therefore do not carry associated antibodies. The exposure of A 0 to leptospires results to animals infected for the first time and still antibody-negative, A 0− , which evolves after a while towards A 2 carrying both IgM and IgG antibodies. Next, A 2 loses IgM antibodies and progresses to A 1 carrying only IgG antibody. When A 1 loses its antibodies, it becomes A 10 again susceptible to secondary and subsequent leptospiral infections. On the other hand, A 1− resulting from an infection of A 1 , evolves towards an A 3 state carrying only IgG but at a high level compared to A 2 and A 1 . And A 3 becomes, A 10 , susceptible again to infections when it loses its antibodies. www.nature.com/scientificreports/ Now, to build the overall model, this antibody kinetics was combined with the infection dynamics according to a SEIRS model as follows. Let N be the constant total number of individuals or animals in the population. The total population is subdivided into ten compartments where (chronologically according to Fig. 1

, right bottom):
• S 0 : number of susceptible animals that have never been neither exposed to nor infected with leptospires. S 0 are antibodies negative. • E 0 : number of exposed animals that have been infected with leptospires for the first time. It is the S 0 that becomes E 0 after exposure and infection with pathogens. E 0 are infected but not infectious yet. E 0 are antibodies negative. • EA 2 : number of animals exposed and infected (not infectious) for the first time and carrying both IgM and IgG antibodies. EA 2 follows after E 0 . • IA 2 : number of animals infected for the first time and now infectious (capable of transmission to susceptible animals) and carrying both IgM and IgG antibodies. IA 2 follows after EA 2 . • IA 1 : number of animals infected for the first time and still infectious but carrying only IgG antibodies. IA 1 follows after IA 2 loses IgM antibodies. • RA 1 : number of animals infected for the first time that have recovered (no longer infectious) and still carrying IgG antibodies. RA 1 follows after IA 1 loses the infection following recovery. • S 1 : number of animals susceptible to leptospiral infection but having already had a history of infection and not carrying associated antibodies. S 1 follows after RA 1 or RA 3 loses the antibodies. • E 1 : number of exposed animals having already had a history of infection and not carrying associated antibodies. It is the S 1 that becomes E 1 after secondary and subsequent infections with pathogens. E 1 are infected but not infectious yet. • IA 3 : number of infected animals (secondary and subsequent infections) that are infectious (capable of transmission to susceptible animals) and carrying high level (greater than in both IA 2 and IA 1 ) IgG antibodies. IA 3 follows after E 1 . • RA 3 : number of infected animals (secondary and subsequent infections) that have recovered (no longer infectious) and still carrying high level IgG antibodies. RA 3 follows after IA 3 loses the infection following recovery.
And, let (t) = (β/N) × [IA 2 (t) + IA 1 (t) + IA 3 (t)] represents the time-dependent force of infection. The combined kinetics of infection and antibodies dynamics, according to Fig. 1, can be described by the set of ten delayed differential equations as:    www.nature.com/scientificreports/ where φ(t k ) = φ k = e −µt k are survival fractions and µ is the natural mortality rate; The summation of all equations equals to zero due to the total population N is kept constant. Each differential equation in Eq. (1) describes the variation over time of the number of animals in the considered compartment. The positive and negative terms after the equal sign in each differential equation represent the increase (incoming arrows in Fig. 1) and the decrease (outgoing arrows in Fig. 1) of the number of individuals, respectively, in the considered compartment.
To keep the total population constant, all mortality (terms " −µ × compartment ") is replaced (term " µN ") by naïve susceptible S 0 . All other terms count for infection and/or transition from one stage to another with transition rates and durations of stays (lag times) in the stages described in Table 1.
The major difference in the dynamics between primo-infected and secondary and later infected manifests itself at two levels: • Primary infected: the class of exposed animals (infected non-infectious) has two populations: those carrying or not antibodies • Secondary and subsequent infected: the class of exposed animals (infected non-infectious) has only one population: that which does not carry antibodies. On the other hand, the level of antibodies in the other classes is higher than for the primary infected.
From the combined models, the meaning of variables for the kinetics of infection according to a SEIRS epidemiological model are given by using the following equations (See S1. Methods for the reduced model), where S, E, I and R are the number of susceptible, exposed, infectious and recovered animals, respectively. The antibody classes are defined from the Fig. 1 as: The basic reproduction number ( R 0 ) is an important epidemiologic metric used to describe the transmissibility of infectious disease. It provides the expected number of secondary cases in a naïve population generated by an infectious animal throughout the infectious period. For this system (involving SEIRS model), the basic reproduction number is given by, R 0 = prob. to be infected during infectious period × fraction of surviving infectious] × population size . This reads as: Note that when β = 0 , the R 0 = 0 , and when β → ∞ , R 0 = αN/(α + µ) . In the absence of any information on β , the R 0 can be inverted to express the contact-transmission rate β as a function of R 0 , population size N and other parameters of the system as:  Table 1, retrieved from the literature (Table S1). Time lags in Eq. (1) for SEIRS (Eq. S1) and antibody kinetics are constrained as shown in Fig. 1 by the following relations: To account for the variability and uncertainty in durations d = t ea , t i , t r and 1/ω 2 in Table 1, the d were all sampled using a Weibull distribution, , with the shape parameter k = 2 and scale parameter, d = 2 √ π × t ea , t i , t r , 1/ω 2 , , where · · · designates the corresponding mean value given in Table 1. The other parameters were calculated using the relations in Eq. (6). Data analysis. The diagnosis of leptospirosis is complicated and depends on the laboratory test. To detect antibodies, there are tests such as microscopic agglutination test (MAT) and immuno-enzymatic assays (ELISA) that are well-known methods. Here, we consider a diagnostic test that would provide one, two or all three of the following outcomes or test results (TR) as shown in Fig. 1: TR 1 -detecting only IgM (and assuming that IgG are also present) and corresponding to the total number of EA 2 and IA 2 as A 2 ; TR 2 -detecting low levels of IgG corresponding to the number of IA 1 and RA 1 as A 1 ; and TR 3 -detecting high levels of IgG corresponding to the number of IA 3 and RA 3 as A 3 .
As shown in Fig. 6, such a diagnostic test provides information regarding the prevalence p i or the fraction of (real or true) positive animals carrying targeted antibodies, the fraction q i of positive animals that are infected, and the lag time ( L 1 and L 2 ) since infection. To extract that information, the reasoning below goes as follows.
• Prevalence of antibody positive: any diagnostic test is characterized by a sensitivity Se and specificity Sp .
Therefore, the fraction of positive T +,i from the diagnostic test result is given by, where p i is the prevalence or the fraction of real or true positive from test i . Now, inverting Eq. (7) provides p i as, www.nature.com/scientificreports/ Then, using outputs of simulations, we can use the determined prevalence p i to infer the value of R 0 . • Fraction of infected: using simulation outputs, we can determine the fraction of infected q i corresponding to the determined R 0 . • Time lag since infection (for a test performed at time t): TR 1 identifies freshly infected animals A 2 , and A 1 and A 3 positive animals should be observed later at L 1 and L 2 , respectively. TR 2 reflects animals A 1 infected about L 1 ago that is expected to be A 3 positive later at L 2 − L 1 . TR 3 identifies animals A 3 infected at least about L 2 ago.
At the steady state, the expressions of p i are given in Eq. (S5). However, we found that simulations of prevalence of antibody positive can be fitted as well with the power law, where a i and k i are given in Table 2.
Likewise, the expressions of the fraction of infected q i at the steady state are given in Eq. (S6) as, Clearly, the q i,s are constant as a function of R 0 .

Simulations details and Statistical analyses.
Stochastic simulations numerically solve Eq. (1) with parameters in Table 1 at daily time step and with initial 10,000 individuals over a 30 year period. Unless stated otherwise, all the parameters are kept the same in all simulations except for the contact-transmission β value that is varied via R 0 using Eq. (5). 1000 simulations were used for statistical analysis. The average numbers of each compartment in Eq. (1) were recorded every 30 days to generate a monthly data series. All simulations were performed using the MATLAB R2016b and statistical analyses were performed using R version 4.0.2.

Data availability
The datasets generated during analyzed and/or the current study were made available from the corresponding author based on reasonable requests.